REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis 
Highway,  Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a 
collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY) 
31-05-2006 

2.  REPORT  TYPE 

Technical  Paper 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

Multi-Domain  Plasma  Expansion  Simulations  Using  a  Particle-in-Cell  Method 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

Lubos  Brieda  (ERC);  Douglas  VanGilder  (Combustion  Research  and  Flow  Technology) 

5d.  PROJECT  NUMBER 

48470052 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Faboratory  (AFMC) 

AFRF/PRSS 

1  Ara  Drive 

Edwards  AFB  CA  93524-7013 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFRL-PR-ED-TP-2006- 1 67 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR’S 
ACRONYM(S) 

Air  Force  Research  Faboratory  (AFMC) 

AFRF/PRS 

5  Pollux  Drive 

Edwards  AFB  CA  93524-70448 

11.  SPONSOR/MONITOR’S 
NUMBER(S) 

AFRL-PR-ED-TP-2006- 1 67 

12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited  (AFRF-ERS-PAS-2006-1 16) 

13.  SUPPLEMENTARY  NOTES 

Presented  at  the  42nd  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference,  Sacramento,  CA,  9-12  July  2006. 

14.  ABSTRACT 

The  statistical  nature  of  the  Particle-In-Cell  (PIC)  algorithm  for  plasma  modeling  requires  that  a  large  number  of  computational 
particles  is  used  per  cell  to  reduce  the  numerical  noise.  This  requirement  presents  a  computational  obstacle  in  cases  involving  rapidly  decaying 
plasmas,  such  as  in  simulations  of  plume  expansion  from  electric  propulsion  (EP)  thrusters.  The  simulation  domain  typically  contains  plasma 
densities  ranging  from  1017  to  1010  particles/m3.  Several  approaches  for  retaining  a  sufficient  per-cell  particle  count  exist,  including  growth  of 
simulation  cells,  particle  splitting,  and  particle  tracking  limited  to  the  back  flowing  particles,  but  none  of  these  is  without  associated  problems. 

In  this  paper,  we  present  an  alternative  approach  based  on  a  multi-domain  modeling.  A  coarse  simulation  is  used  to  sample  particle  flux  into  a 
subdomain  enclosing  the  region  of  interest.  A  second  simulation  is  then  performed  on  the  subdomain,  with  particles  injected  at  domain 
boundaries  according  to  the  prescribed  flux.  This  approach  is  used  to  predict  ion  current  to  a  simple  cylindrical  probe  located  on  a  satellite  using 
a  cluster  of  four  Hall  thrusters  for  primary  propulsion.  The  effect  of  sheath  resolution  is  investigated  and  results  are  compared  to  an  analytical 
model. 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE 

OF  ABSTRACT 

OF  PAGES 

PERSON 

Daron  Bromaghim 

a.  REPORT 

b.  ABSTRACT 

c.  THIS  PAGE 

A 

13 

19b.  TELEPHONE  NUMBER 

(include  area  code) 

Unclassified 

Unclassified 

Unclassified 

N/A 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std.  239.18 


Unclassified 


42nd  AIAA/ASME/SAE/ASEE  Joint  Propulsion  Conference  &  Exhibit 
9-12  July  2006,  Sacramento,  California 


AIAA  2006-5023 


Multi-Domain  Plasma  Expansion  Simulations  Using  a 

Particle-in-Cell  Method 

Lubos  Brieda* 

Air  Force  Research  Laboratory,  Edwards  AFB,  CA  93524 

and 

Douglas  VanGilder* 1' 

Combustion  Research  and  Flow  Technology,  Pipersville,  PA  18947 


The  statistical  nature  of  the  Particle-In-Cell  (PIC)  algorithm  for  plasma  modeling  re¬ 
quires  that  a  large  number  of  computational  particles  is  used  per  cell  to  reduce  the  numer¬ 
ical  noise.  This  requirement  presents  a  computational  obstacle  in  cases  involving  rapidly 
decaying  plasmas,  such  as  in  simulations  of  plume  expansion  from  electric  propulsion  (EP) 
thrusters.  The  simulation  domain  typically  contains  plasma  densities  ranging  from  1017  to 
1010  particles/m3.  Several  approaches  for  retaining  a  sufficient  per-cell  particle  count  exist, 
including  growth  of  simulation  cells,  particle  splitting,  and  particle  tracking  limited  to  the 
backflowing  particles,  but  none  of  these  is  without  associated  problems.  In  this  paper,  we 
present  an  alternative  approach  based  on  a  multi-domain  modeling.  A  coarse  simulation 
is  used  to  sample  particle  flux  into  a  subdomain  enclosing  the  region  of  interest.  A  second 
simulation  is  then  performed  on  the  subdomain,  with  particles  injected  at  domain  bound¬ 
aries  according  to  the  prescribed  flux.  This  approach  is  used  to  predict  ion  current  to  a 
simple  cylindrical  probe  located  on  a  satellite  using  a  cluster  of  four  Hall  thrusters  for  pri¬ 
mary  propulsion.  The  effect  of  sheath  resolution  is  investigated  and  results  are  compared 
to  an  analytical  model. 


I.  Introduction 

The  Air  Force  Research  Laboratory  (AFRL)  at  Edwards  AFB  is  developing  a  computational  package 
called  COLISEUM1  for  modeling  the  interaction  of  electric  propulsion  (EP)  plumes  with  spacecraft  compo¬ 
nents.  COLISEUM  consists  of  three  primary  modules:  RAY,  AQUILA  and  DRACO.  The  RAY  module  is 
used  to  quickly  compute  line-of-sight  contamination  predictions,  while  AQUILA2  and  DRACO3  are  particle 
codes  that  track  plume  expansion  from  the  source  using  the  Particle-In- Cell  (PIC)4  method  combined  with 
MCC/DSMC5,6  treatment  of  collisions. 

A.  Density  Decay  in  Particle  Codes 

The  kinetic  nature  of  the  PIC  and  MCC/DSMC  methods  requires  that  a  large  number  of  particles  is  used 
per  computational  cell  to  reduce  the  statistical  noise  and  to  resolve  the  local  velocity  distribution  function 
(VDF).  Due  to  computational  limitations,  plasma  is  described  by  macroparticles,  with  each  macroparticle 
representing  w  real  ions.  As  shown  in  this  paper  and  summarized  in  table [I]  the  plasma  density  around  a  typ¬ 
ical  EP  spacecraft  in  GEO  ranges  from  1017  particles/m3  near  the  thruster  exit  to  around  1010  particles/m3 
in  shadowed  areas.  Retaining  a  sufficient  particle  count  per  cell,  N ,  over  a  density  decay  of  seven  orders  of 
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magnitude  is  problematic.  If  the  simulation  cell  size  remains  constant,  N  scales  with  density  as 


n2  = 

m 

region 

n(#/m-3) 

HET  exit 

1017 

CEX  wing 

1012  -  1014 

Wake 

1010 

Table  1.  Plasma  densities  encountered  around  a  typical  spacecraft  flying  in  GEO  and  using  a  Hall  Thruster 
for  propulsion. 

Retaining  7Vm^n  =  2  (minimum  for  a  DSMC  collision)  in  the  wake  region,  requires  the  simulation  to  track 
over  10  million  particles  per  cell  in  the  high  density  region  near  the  thruster  exit.  Even  if  the  wake  region  is 
deemed  insignificant,  and  the  simulation  is  limited  to  a  study  of  CEX  plasma  interacting  with  solar  arrays, 
the  simulation  needs  to  track  around  100,000  particles  per  cell  in  the  high  density  region.  The  requirement 
of  two  particles  per  cell  is  not  sufficient  to  properly  resolve  the  VDF,  and  could  lead  to  erroneous  predictions 
of  surface  sputtering. 

B.  Approach  1:  Cell  Growth 

A  commonly  used  approach  in  treating  the  density  decay  is  allowing  the  cell  size  to  grow  inversely  with  the 
decay  in  density.  This  approach  is  used  by  AQUILA.  While  it  achieves  a  more  uniform  per-cell  macroparticle 
counts,  the  cell  growth  limits  the  resolution  of  the  solution  in  the  low  density  region.  This  is  especially 
problematic  if  a  detailed  map  of  surface  erosion  is  desired.  Furthermore,  unstructured  codes  typically  use 
a  first-order  potential  solver  (Boltzmann-inversion  or  linear  Poisson  solver),  and  differentiation  of  the  basis 
function  results  in  a  constant  electric  field  per  cell.  This  discontinuous  nature  of  the  field  solution  then 
results  in  an  increased  ion  diffusion  to  walls  due  to  a  non-physical  sheath  widening. 

C.  Approach  2:  Particle  Splitting 

An  alternative  to  the  use  of  an  unstructured  mesh  is  particle  splitting.  In  this  approach,  macroparticles 
entering  a  low  density  region  are  split  into  two  sub-particles  with  a  halved  specific  weight.  Although  this 
approach  increases  the  number  of  particles  per  cell,  it  does  not  necessarily  improve  the  representation  of  the 
local  VDF.  The  collision  cross-section  for  the  CEX  collision  in  Xenon  is  given  by  Pullins[7]  as 

acex  ( Xe,Xe+ )  =  1.1872  x  1CT20  [23.3  log^)  +  188.81]  (2) 

and  collision  mean  free  path  scales  as 

^ m  1/^n^"  (3) 

The  magnitude  of  relative  velocity  in  the  charge  exchange  wing  can  be  approximated  from  neutral  temper¬ 
ature  near  the  thruster  exit.  Assuming  the  neutrals  leave  the  thruster  at  1000K  O.leV),  the  collision 
cross-section  is  1.5  x  10-18,  leading  to  Am  ~  106m.  Hence,  the  split  particles  will  travel  through  the  low 
density  region  on  top  of  each  other,  which  is  analogous  to  the  motion  of  the  original  macroparticle.  Small 
variation  in  velocity  can  be  achieved  by  offsetting  the  new  born  particles  from  the  original  center-of-mass 
trajectory.  However,  the  effect  of  such  random  particle  jumps  on  the  simulation  results  in  not  quantified. 
Alternatively,  the  velocity  of  the  new  born  particles  could  be  adjusted  such  that 

m  (v[  +  v'2)  =  mv  (4) 

The  difficulty  here  lies  in  choosing  the  new  velocities.  If  the  collective  velocity  of  the  incoming  macroparticles 
can  be  represented  by  a  double-Maxwellian  function,  then  the  two  velocity  components  could  be  chosen 
from  the  drift  and  thermal  components  of  /m,i  and  /m,2-  This  approach  requires  each  particle  to  track  a 
representative  VDF,  which  will  vary  as  a  function  of  position.  Such  a  detailed  treatment  of  particle  motion 
is  beyond  the  scope  of  standard  PIC  codes. 
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D.  Approach  3:  Fluid  Beam  Model 

Additional  alternative  approaches  were  explored,  of  these  probably  the  most  notable  is  the  division  of  the 
plume  into  the  primary  ion  beam  and  the  CEX  wing  region,  as  done  by  Roy [8]  and  Wang [9].  The  main 
beam  is  described  by  an  analytical  profile.  The  kinetic  treatment  is  limited  to  the  charge  exchange  ions,  and 
a  lower  macroparticle  weight  can  be  used.  This  approach  improves  the  velocity  representation  in  the  low 
density  region,  however,  it  is  not  mass-conserving.  This  simplification  could  lead  to  an  artificially  high  ion 
currents  in  vacuum  chamber  simulations,  in  which  loss  of  fast  ions  to  CEX  collisions  is  not  negligible  due 
to  the  increased  background  neutral  pressure.  Furthermore,  the  main  beam  is  described  by  an  analytical 
profile,  which  is  at  best  an  approximation  based  on  a  witnessed  beam  divergence. 

II.  Multi-Domain  PIC  Model  (MD-PIC) 


A.  Formulation 

This  paper  presents  a  new  approach  for  resolving  the  ion  dynamics  in  a  low-density  region  based  on  a  multi- 
domain  formulation.  The  simulation  is  performed  in  three  steps,  with  the  first  two  steps  used  to  collect  flux 
of  ions  into  a  region  of  interest.  Then,  a  detailed  simulation  is  performed  on  the  subdomain.  This  approach 
is  applicable  in  cases  where  fine  details  of  plasma  dynamics  in  the  subdomain  region  do  not  have  a  strong 
coupling  to  the  bulk  plasma  expansion.  The  three  steps  used  in  the  process  are  described  below. 

•  Plume  Expansion:  The  simulation  first  computes  the  plume  expansion  from  the  thruster  using  a 
coarse  mesh  and  a  simplified  geometry  description.  This  step  is  identical  to  a  typical  plume/spacecraft 
interaction  modeling.  The  quasi- neutral (QN)  assumption  generally  holds,  and  the  potential  can  be 
computed  from  the  Boltzmann  inversion, 


4>  =  4>o  +  kTe  In  (rii/no) 


(5) 


Collisions  are  treated  using  the  MCC  or  DSMC  approach.  The  simulation  continues  until  a  steady 
state  is  achieved,  dM/dt  ~  0  and  dp/d t  ~  0,  where  M  is  the  number  of  macroparticles  and 


M 

^ rnjVi 


(6) 


is  the  average  particle  momentum. 

•  Flux  Sampling:  Next,  particle  flux  into  the  region  of  interest  is  sampled  over  a  large  number  of  time 
steps.  The  region  of  interest  is  typically  located  in  a  low  density  area,  and  the  long  sampling  time  is 
necessary  to  accurately  capture  the  VDF  of  the  incoming  particles.  A  discrete  approach  is  used  in  this 
method.  Instead  of  attempting  to  qualify  the  incoming  flux  using  an  analytical  profile  (summation  of 
Maxwellian  functions,  etc. . . ),  the  code  simply  outputs  the  final  location  and  velocity  of  each  incoming 
particle  into  a  file.  Mass  flow  rate  of  the  sampled  particles  is  computed  from 

mMsw 


where  Ms  is  the  number  of  sampled  particles. 

•  Detailed  Modeling  The  last  step  involves  detailed  modeling  of  the  low-density  region.  The  simulation 
domain  is  limited  to  the  subdomain  enclosing  the  region  of  interest.  Particles  are  injected  into  the 
domain  using  a  volume  source  which  returns  a  random  particle  from  the  list  sampled  in  the  previous 
step.  Number  of  particles  to  sample  each  time  step  is  computed  from  m  calculated  in  the  previous  step. 
Limiting  the  simulation  to  the  small  subdomain  provides  a  two-fold  benefit.  First,  since  the  simulation 
describes  a  smaller  physical  volume,  a  finer  cell  size  can  be  used  while  retaining  a  computationally 
feasible  number  of  mesh  nodes.  The  reduction  in  cell  size  results  in  an  increased  spatial  resolution,  and 
hence  a  more  detailed  surface  mesh  can  be  used.  Second,  the  simulation  tracks  only  the  particles  in  the 
subdomain.  The  macroparticle  specific  weight  can  thus  be  reduced  by  several  orders  of  magnitudes, 
resulting  in  a  high  per-cell  macroparticle  count. 
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B.  Simulation  Setup 


The  MD-PIC  approach  is  demonstrated  in  this  paper  by  modeling  ion  current  collection  by  a  Langmuir 
probe  located  on  a  hypothetical  GEO  spacecraft.  The  approach  is  used  to  investigate  the  importance  of 
resolving  the  plasma  sheath  on  current  predictions.  The  spacecraft  and  the  simulation  domain  are  shown  in 
figure [l(a)|  The  spacecraft  uses  a  cluster  of  four  200W  Hall  Effect  Thrusters10  for  propulsion.  Each  thruster 
is  operating  at  0.8 A  discharge  current  and  250V  potential.  A  1.5m  long  beam  extends  from  one  side  of 
the  spacecraft.  At  the  end  of  the  beam  is  located  a  10cm- long  cylindrical  probe.  Due  to  limits  of  mesh 
resolution,  the  diameter  of  the  probe  was  set  to  5cm  in  the  coarse  simulation,  which  matches  the  diameter 
of  the  holding  beam.  The  probe  diameter  was  decreased  to  1.2cm  in  the  subdomain  model.  The  extent  of 
the  subdomain  is  outlined  by  the  purple  boundary.  Close  up  of  the  beam  tip  and  the  attached  probe  on  the 
fine  mesh  is  shown  in  figure [1(b)]  This  plot  also  illustrates  the  ratio  in  mesh  size  between  the  coarse  and  the 
fine  mesh. 


(a)  Simulation  Domain  (b)  Sub-Domain  Mesh  Detail 


Figure  1.  Simulation  Setup.  The  hypothetical  satellite  contains  an  array  of  four  HET  thrusters,  and  a  long 
beam  holding  a  lOcm-long  Langmuir  probe.  The  region  enclosed  by  the  sub-domain  is  shown  in  purple.  Figure 
(b)  shows  a  close  up  of  the  probe  on  the  fine  mesh. 


1.  Plume  Expansion 


Step  1  (plume  expansion)  was  modeled  using  a  uniform  70  x  123  x  100  mesh  with  Ax  =  0.05m.  Symmetry 
was  assumed  along  the  z  —  y  plane,  and  only  half-domain  was  simulated.  A  reflective  particle  boundary 
condition  was  enforced  on  the  plane  of  symmetry.  The  remaining  faces  were  open  and  particles  leaving 
through  these  faces  were  removed  from  the  simulation.  Potential  was  computed  using  the  quasi-neutral 
approximation  (eq.  [5|,  with  a  constant  kTe  =  1.5eV.  Conditions  at  thruster  exit  were  used  to  set  <po  =  20V 
and  no  =  8  x  1016m-3.  The  spacecraft  potential  was  set  to  0V,  but  the  thruster  was  assumed  to  float  at 
-f  20V.  The  simulation  ran  for  10,000  time  steps  with  At  =  5  x  10-7s.  Total  of  3.8  million  particles  were 
tracked  at  the  end  of  the  simulation. 

Ions  were  injected  using  a  particle  source  based  on  LIF11  velocity  distribution.  The  specific  weight  of 
ions  was  5  x  108.  Neutral  mass  flow  rate  was  approximated  as  rhn  =  O.lra*.  The  simulation  did  not  track  the 
neutral  particles,  instead  a  steady-state  distribution  was  modeled  by  loading  a  background  neutral  density 
by  projecting  the  source, 

nn  =  w  (8) 


where 


/ 


vfM  d3u 


(9) 
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The  CEX  collisions  were  treated  using  the  MCC  approach.  For  every  particle,  collision  frequency  was 
calculated  as 

v  =  nnacExg  (10) 

with  g  =  Vi,  implying  that  neutral  velocity  is  insignificant  compared  to  the  ion  velocity.  Collision  probability, 

P  =  1  —  exp  (— vAt)  (11) 

was  then  compared  to  a  random  number,  R,  and  the  collision  was  performed  for  R  <  P.  The  charge  transfer 
was  modeled  by  replacing  the  ion  velocity  with  a  random  neutral  thermal  velocity, 

j2kTn 

vth  =  \  - /m  (12) 

V  m 

where  Tn  =  500 K  is  a  constant  background  neutral  temperature. 

2.  Sampling 

The  particle  flux  into  the  subdomain  was  collected  for  additional  50,000  time  steps.  The  dimensions  of  the 
subdomain  were  0.4m x  0.4m x  0.6m  and  the  subdomain  was  centered  in  the  x  —  z  plane  on  the  probe.  The 
sampling  resulted  in  collection  of  74893  macroparticles,  at  mm  3.267  x  10  10kg/s. 

3.  Detailed  Simulation 

The  collected  particle  flux  was  then  used  to  inject  particles  into  the  subdomain.  This  simulation  used  a  fine 
mesh  with  Ax  =  0.005m  and  80  x  80  x  120  cells.  A  detailed  representation  of  the  probe  with  d  =  1.2cm  was 
used  in  this  case.  The  specific  weight  was  decreased  by  0(4)  to  5  x  104  and  the  subdomain  contained  around 
2  million  particles  at  steady  state.  The  effect  of  sheath  resolution  on  simulated  collected  current  was  studied 
by  computing  the  potential  using  two  methods:  inversion  of  the  Boltzmann  equation  (QN  approximation) 
and  by  solving  the  Poisson’s  equation.  An  I  —  V  curve  was  generated  for  both  potential  solvers  by  biasing 
the  probe  to  six  potential  values,  fo:  (j)p  —  30V,  <j)p  —  20V,  <j)p  —  10V,  <fip  +  10V  and  rf>p  +  20V,  where 
c j)p  =  3V  is  the  simulation  plasma  potential  in  the  vicinity  of  the  probe.  Each  simulation  was  executed  for 
5000  time  steps  with  At  =  1  x  10-7m/s. 


III.  Results 


Results  on  the  Coarse  Mesh 


Figure  [2]  shows  the  potential  profile  around  the  spacecraft  obtained  after  5000  time  steps  on  the  coarse  mesh. 
The  potential  ranges  from  20V  at  the  thruster  exit  to  -24V  in  plasma-free  region.  The  potential  near  the 
probe  tip  is  3V.  Despite  the  simulation  tracking  nearly  4  million  particles,  the  potential  profile  is  very  noisy, 
especially  in  the  wake  on  the  front  side  of  the  spacecraft,  and  around  the  probe.  The  noise  is  directly  related 
to  the  density  decay  and  corresponding  low  macroparticle  count,  which  are  plotted  in  figures [3(a)]  and [3(b) [ 
As  predicted  in  the  introduction,  number  of  macroparticles  per  cell  decreased  from  20,000  near  the  thruster 
to  below  1  at  the  far  extent  of  the  CEX  wing  and  around  the  probe  tip. 


As  can  be  seen  from  figure  4(a) ,  these  results  correspond  to  the  steady-state  solution,  however,  the  low 
macroparticle  count  makes  it  very  difficult  to  use  the  simulation  results  to  accurately  predict  the  ion  current 
collection  by  the  probe.  The  ion  current  was  computed  by  dividing  the  collected  current  by  the  probe  area, 


J  = 


I 


2irrl  +  7rr2 


(13) 


where  r  is  the  probe  radius  and  l  is  0.1m.  An  attempt  was  made  to  smooth  the  collected  data,  however,  as 
can  be  seen  from  the  purple  markers,  smoothing  did  not  produce  any  improvement  in  the  data. 
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(a)  side  view,  y  —  z  plane 


(b)  top  view,  x  —  z  plane 


Figure  2.  Potential  profile  on  the  coarse  mesh.  The  simulation  noise  due  to  insufficient  number  of  particles 
per  cell  is  clearly  visible. 


(a)  Ion  Density 


(b)  Particles  per  cell 


Figure  3.  Top  view  of  the  x  —  z  plane  showing  the  ion  density  and  per-cell  particle  count.  Rapid  decay  in 
density  results  in  the  simulation  tracking  less  than  1  particle  per  cell  near  the  end  of  the  beam  holding  the 
Langmuir  probe. 
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(b)  Current  Collection 


Figure  4.  Figure  (a)  shows  the  progression  of  particle  count  and  average  momentum  versus  time  step.  Both 
properties  reach  a  constant  value  around  it  =  3000,  indicating  a  steady-state.  Figure  (b)  shows  the  current 
collection  by  the  probe  versus  time  step.  The  results  are  very  noisy  due  to  the  low  particle  count,  and 
smoothing,  shown  by  the  purple  markers,  does  not  improve  the  results. 


B.  MD-PIC  Results 

The  continuity  of  solution  was  inspected  visually,  as  shown  in  figure  [5]  The  left  plots  show  the  potential 
and  ion  density  obtained  on  the  coarse  mesh.  The  outline  of  the  region  covered  by  the  subdomain  is  also 
outlined.  The  right  plots  show  the  corresponding  solution  on  the  subdomain.  The  plots  show  instantaneous 
values,  obtained  at  corresponding  steady-states,  which  lead  to  some  variation  in  values  around  the  border. 
However,  both  solutions  display  the  same  trends,  indicating  that  continuity  is  preserved. 

The  coarse  results  indicate  a  significant  amount  of  numerical  noise.  The  noise  is  due  to  both  the  large 
simulation  cell  size,  and  the  low  per-cell  macroparticle  count.  A  drastic  improvement  can  be  noticed  in  the 
MD-PIC  results.  These  results  show  a  clear  wake  behind  the  probe,  feature  which  is  not  distinguishable  in 
the  coarse  simulation.  The  density  wake  results  in  a  potential  gradient  around  the  probe.  Although  a  similar 
drop  can  be  noticed  in  the  coarse  result,  it  is  not  clear  whether  the  drop  is  a  steady-state  prediction,  or  just 
a  transient  feature  existing  due  to  numerical  noise.  The  gradient  also  extends  over  a  larger  physical  area, 
which  is  a  direct  consequence  of  the  excessively  large  simulation  cells. 

C.  Effect  of  Sheath  Resolution 

The  quasi-neutral  approximation  directly  inverts  the  density  at  each  node  without  any  consideration  to 
potential  on  surrounding  nodes.  As  such,  it  screens  out  any  imposed  potential,  and  hence  artificially  limits 
the  size  of  a  sheath  to  the  length  of  the  simulation  cell  size.  The  DADI  method,  on  the  other  hand,  solves 
the  Poisson’s  equation, 

V2</>  =  -  —  (14) 

£o 

and  does  resolve  the  sheath.  The  disadvantage  of  the  DADI  method  is  the  additional  amount  of  time  needed 
to  converge  the  iterative  solver.  The  effect  of  the  sheath  resolution  on  current  collection  was  studied  by 
comparing  the  solutions  obtained  using  the  QN  approximation  to  the  DADI  solution.  Comparison  of  plasma 
potentials  and  densities  is  shown  in  strips  |6(a)|6(f)j  and  |7(a)]|7(f)|  The  top  plots  show  the  results  obtained 
using  the  QN  solver,  while  the  bottom  half  contains  the  results  from  the  DADI  runs.  Both  solvers  result  in 
similar  potential  profiles,  but  the  DADI  solutions  are  smoother.  The  smoothness  is  due  to  the  elliptic  nature 
of  the  Poisson’s  equation.  Important  to  notice  is  the  lack  of  a  sheath  around  the  probe  in  the  QN  results  for 
(f)b  =  (ftp  —  20H  and  fa  =  fa  +  20V. 

The  failure  to  resolve  the  sheath  leads  to  noticeable  differences  in  density  profiles.  The  QN  approach 
results  in  densities  nearly  invariant  with  the  applied  probe  potential,  and  wake  formation  is  primarily  due 
to  the  physical  obstruction  of  the  flow.  On  the  other  hand,  the  wake  structure  varies  greatly  in  the  DADI 
results.  The  potential  applied  to  the  probe  extends  to  the  plasma  through  the  sheath,  and  increases  the 
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(a)  Potential,  Coarse  Mesh 


(c)  Ion  Density,  Coarse  Mesh 


(b)  Potential,  Subdomain 


(d)  Ion  Density,  Subdomain 


Figure  5.  Continuity  of  solution.  The  plots  show  the  instantaneous  values  at  corresponding  steady-states, 
which  lead  to  some  variation  along  the  domain  boundary. 
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(d)  DADI,  cf)b  =  (fip  —  20V 


(e)  DADI,  cf)b  =  (f)p  +  0V  (f)  DADI,  cf)b  =  0P  +  20V 


Figure  6.  Potential  on  the  subdomain,  obtained  using  the  QN  and  DADI  methods,  for  several  probe  bias 
potentials. 
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Figure  7.  Ion  density  on  the  subdomain,  obtained  using  the  QN  and  DADI  methods,  for  several  probe  bias 
potentials. 
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effective  probe  area.  Large  negative  bias  leads  to  an  increase  in  plasma  density  around  the  probe.  The  wake 
is  almost  non-existent,  as  the  slow  moving  CEX  ions  are  attracted  back  to  the  probe.  Important  to  note  is 
the  similarity  between  the  two  solvers  for  the  fa  =  <j>p  case.  This  result  was  expected,  since  at  this  condition 
the  probe  is  biased  to  the  plasma  potential,  and  subsequently,  there  is  no  sheath.  The  positive  probe  bias 
results  in  repulsion  of  ions  from  the  probe.  The  probe  is  biased  to  -f  20V  from  the  surrounding  plasma,  and 
hence  only  particles  with  KE>20V  have  sufficient  energy  to  reach  it.  Since  the  CEX  ions  are  born  nearly 
at  rest  (neutral  thermal  velocity),  and  the  potential  drop  between  the  thruster  and  the  probe  is  <17V,  this 
bias  should  result  in  no  current  collection. 

This  is  indeed  the  case,  as  can  be  seen  from  the  ion  current  plots  in  figure [8(a)]  The  green  line  indicates  the 
current  collected  versus  the  time  step  for  fa  =  <j)p  +  20V.  No  statistically  significant  current  is  sampled.  As 
expected,  an  overlap  is  seen  between  the  black  dotted  and  the  black  solid  lines.  The  dotted  lines  correspond 
to  the  DADI  results,  while  the  solid  line  shows  the  current  obtained  from  the  QN  approximation.  The  black 
lines  indicate  the  ion  current  at  the  floating  potential,  fa  =  fa.  Since  this  condition  does  not  result  in 
formation  of  a  sheath,  the  lack  of  sheath  resolution  by  the  QN  method  is  trivial.  A  strong  dependence  on 
sheath  resolution  is  seen  in  the  plot  for  the  negative  probe  bias.  The  potential  gradient  in  the  sheath  modifies 
the  trajectory  of  particles  that  would  otherwise  flow  past  the  probe,  and  hence  the  DADI  case  results  in  a 
higher  current  collection.  The  results  are  summarized  in  the  I  —  V  curve  plotted  in  figure  |8(b)|  As  can  be 
seen,  the  probe  starts  repelling  ions  at  fa  >  10V.  The  QN  case  under-predicts  the  current  due  to  screening 
of  applied  potential  at  length  scales  larger  than  the  cell. 

The  two  MD-PIC  results  were  also  compared  to  a  theoretical  model  of  Prokopenko  and  Laframboise,12 
given  by 


jpl  =  jo 


2  (Q/7r)1/2  +  exp  (Q)  erfc  {Q1/2^j 


(15) 


where 

Q  =  ~qfa/kTe  (16) 


This  relationship,  derived  assuming  an  infinitely  long  cylinder,  the  thick- sheath  limit  and  uniform  mono- 
energetic  beam  at  infinity,  scales  with  kTe  and  jo.  The  comparison  was  performed  using  kTe  =  1.5eV 
and  jo  =  Jqjv(+0V).  An  agreement  within  an  order  of  magnitude  is  seen  between  the  model  and  both 
simulation  results,  however,  the  DADI  case  results  in  a  significant  error  reduction,  and  a  comparable  I  vs. 
V  behavior.  The  trend  seen  in  the  QN  curve  corresponds  to  kTe  of  12eV.  The  disrepancies  in  magnitude 
could  be  attributed  to  a  limited  applicability  of  this  model  to  the  simulated  case. 


Figure  8.  Probe  current  collection  determined  by  the  MD-PIC  approach.  Figure  (b)  shows  the  I  —  V  curve 
for  the  ion  current. 


IV.  Conclusion  and  Future  Work 

A  multi-domain  approach  (MD-PIC)  for  resolving  ion  dynamics  in  low-density  regions  was  presented  in 
this  paper.  This  approach  is  based  on  a  three  step  modeling  approach,  with  the  first  two  steps  used  to  sample 
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particle  flux  into  a  region  of  interest.  Detailed  simulation  on  the  subdomain  is  then  performed,  utilizing  a 
reduced  cell  size,  refined  surface  geometry,  and  a  lower  macroparticle  specific  weight.  The  MD-PIC  approach 
was  used  in  this  paper  to  study  the  effect  of  sheath  resolution  on  collection  of  ion  current  by  a  Langmuir 
probe  located  on  a  hypothetical  EP  GEO  spacecraft.  An  I  —  V  curve  was  obtained  using  six  probe  bias 
potentials.  Current  collection  demonstrated  a  strong  dependence  on  the  sheath  resolution  for  all  negative 
bias  potentials.  The  quasi-neutral  case,  which  does  not  resolve  the  sheath,  under-predicted  the  current 
collection,  due  to  the  reduced  collection  area  associated  with  the  numerical  screening  of  potential  gradient 
at  length  scales  larger  than  the  cell.  Simulation  results  were  also  compared  to  a  theoretical  model,  and  an 
agreement  within  30%  was  achieved  for  the  Poisson  solver  case. 

These  results  indicate  that  sheath  resolution  is  critical  in  cases  where  a  detail  prediction  of  surface  ion 
collection  is  needed.  The  sheath  increases  the  effective  collection  area,  and  hence  increases  the  collected 
particle  flux.  The  potential  profile  within  the  sheath  also  alters  the  trajectory  of  the  incoming  particles. 
Correct  prediction  of  incident  angle  is  important  in  sputtering  analysis,  as  sputter  yield  shows  a  strong 
dependence  on  incident  angle.  The  MD-PIC  approach  provides  the  capability  to  resolve  the  sheath  in  the 
region  of  interest,  while  retaining  a  high  per-cell  particle  count.  As  such,  this  approach  is  applicable  to 
studies  of  surface  erosion,  spacecraft  charging,  and  detailed  modeling  of  plasma  probes.  In  the  future,  an 
attempt  will  be  made  to  further  validate  the  model  through  a  comparison  with  published  in-flight  sensor 
measurements. 
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